rm(list=ls())
library(logr)
library(ggplot2)
#setwd("/Users/tg2778/Dropbox/0_Reviews_RnRs/072022_JOP_Roads/v3_JOP/Replication - Roads")
log_open("../4_Log/01_RDpolotLog.log")
load("../1_Data/censusvillageacroads.Rdata")

censusvillageacroads$pmgsyroad<-0
censusvillageacroads$pmgsyroad[censusvillageacroads$roadcode>0 & is.na(censusvillageacroads$roadcode)==FALSE]<-1

censusvillageacroads$popcutoff500<-censusvillageacroads$pc01_vd_t_p - 500

census84<-subset(censusvillageacroads,censusvillageacroads$popcutoff500>=-86 & censusvillageacroads$popcutoff500<=86)

ggplot(census84, aes(popcutoff500, pmgsyroad, popcutoff500 > 0)) +
  stat_summary_bin(aes(group = popcutoff500 > 0), geom = 'point', fun = 'mean') +
  geom_smooth(aes(group = popcutoff500 > 0), se = T, color = 'black', 
              span = 1, method = 'loess') +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  ylab('Probability of PMGSY connectivity by 2018') +
  xlab('Normalized population') +
  theme(strip.text.x = element_text(size = 14))+
  theme(strip.text.y = element_text(size = 14))+
  theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                    panel.background = element_blank(), axis.line = element_line(colour = "black"))
ggsave("../3_Figs/Fig3A.pdf", width=8, height=4.5)
log_close()